package scu.maqiang.piezoelectricity;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES2T3PE;
import scu.maqiang.mesh.Mesh2T3;
import scu.maqiang.mesh.T3Type;
import scu.maqiang.numeric.DirectSSolver;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.SRMatrix;

public class Piezoelectricity2DT3 {

	public static void main(String[] args) {
		long begin = System.nanoTime();
		Mesh2T3 mesh = new Mesh2T3().square2D(10, 10, T3Type.Left);
		FES2T3PE fs = new FES2T3PE(mesh);
		SRMatrix A = new SRMatrix(fs.GetNdof());
		double[] RHS = new double[fs.GetNdof()];
		double[] x = RHS.clone();
		double[][] matD = new double[][] {{14.9e10, 9.8e10,  0.0,    0.0,     -5.4}, 
			                              {9.8e10,  14.3e10, 0.0,    0.0,     -5.4}, 
			                              {0.0,     0.0,     2.2e10, 11.7,    0.0}, 
			                              {0.0,     0.0,     -11.7,  8.0e-9,  0.0}, 
			                              {5.4,    5.4,    0.0,    0.0,     7.2e-9}};
	    fs.assembleStiff(new double[][][] {matD}, BVPType.CONSTITUTE_MATRIX_COMMON, A);
	    fs.assembleSource(new double[] {100.0, 0.0, 200}, BVPType.CONSTITUTE_MATRIX_COMMON, RHS);
	    fs.applyBC_MBN(A, RHS, new int[] {0, 1}, new double[] {0.0, 0.0}, 1, 2, 3, 4);
	    fs.applyBC_MBN(A, RHS, new int[] {2}, new double[] {0.0}, 1, 2, 3);
	    
	    DirectSSolver dsolver = new DirectSSolver(A);
	    //double[] RHS0 = MVO.copy(RHS);
	    //double[] RHS0 = new double[]
		x = dsolver.LUSolve(RHS, 0, 1e-6);
	    //A.save("SRMatrixFile.txt");
	    //MVO.saveArray("RHSFile.txt", RHS);
	    //Arrays.fill(x, 0.0);
//	    NewIterSSolver solver = new NewIterSSolver(A);
//	    //solver.setEps(1.0e-10);
//	    solver.PCGNE_M(RHS, x, 2);
	    //.GMRES_SSOR(RHS, x, 1.5, 200, 2);
	    double[][] uvphi = new double[3][mesh.getNv()];
	    fs.extractComponent(x, uvphi);
	    mesh.toTecplot("PESolution2D.dat", uvphi);
	    System.out.println("SOlution: ");
	    System.out.println("u min:" + MVO.min(uvphi[0]) + "\tmax: " + MVO.max(uvphi[0]));
	    System.out.println("v min:" + MVO.min(uvphi[1]) + "\tmax: " + MVO.max(uvphi[1]));
	    System.out.println("phi min:" + MVO.min(uvphi[2]) + "\tmax: " + MVO.max(uvphi[2]));
		long end = System.nanoTime();
		System.out.println("Solve time: " + (end - begin) / 1.0e9 + " s");
	}

}
